library(readr)
library(dplyr)
library(stargazer)
library(pscl)
library(lmtest)
library(sandwich)
library(MASS)
library(ggplot2)
library(xtable)
library(betareg)
library(purrr)

bwgdata <- read_csv("bwg_data.csv")
bwgdata$region <- as.factor(bwgdata$region)

###############################################
##### Summary Statistics of Key Variables #####
###############################################

vars <- c("mps_count", "criticism", "praise",
          "social_ties", "elec_comp", "execnat", "log_gdp", "log_milex",
          "log_milper", "log_pop", "gwf_military", "mid", "civil_war",
          "coups_cumulative", "purge", "media_diverse", "hhi")


vars_of_interest <- c("mps_count",
                      "criticism", "praise",
                      "social_ties", "elec_comp", "execnat", "log_gdp", "log_milex",
                      "log_milper", "log_pop", "gwf_military", "mid", "civil_war",
                      "coups_cumulative", "purge", "media_diverse", "hhi")

summary_stats <- map_dfr(vars_of_interest, function(var) {
  x <- bwgdata[[var]]
  tibble(
    variable = var,
    mean  = mean(x, na.rm = TRUE),
    median = median(x, na.rm = TRUE),
    sd    = sd(x, na.rm = TRUE),
    var   = var(x, na.rm = TRUE),
    n     = sum(!is.na(x))
  )
})

# For LaTeX #

summary_stats_round <- summary_stats %>%
  mutate(
    mean   = round(mean, 3),
    median = round(median, 3),
    sd     = round(sd, 3),
    var    = round(var, 3)
  )

print(
  xtable(summary_stats_round,
         caption = "Summary Statistics of Key Variables",
         label = "tab:summary_stats"),
  include.rownames = FALSE,
  booktabs = TRUE,
  sanitize.text.function = identity
)

##################
##### Models #####
##################

#########################
### Table 2 - Model 1 ###
#########################

model_zinb <- zeroinfl(
  mps_count ~ social_ties + region + purge + gwf_military +
    log_gdp + log_pop + mid + civil_war + coups_cumulative +
    log_milex + log_milper + hhi |
    mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb)
coeftest(model_zinb, vcov = robust_vcov)

##########################
### Table 2 - Model 2 ####
##########################

model_zinb2 <- zeroinfl(
  mps_count ~ social_ties * elec_comp + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi | mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov2 <- sandwich::vcovOPG(model_zinb2)
coeftest(model_zinb2, vcov = robust_vcov2)

# Prediction #

pred_data <- expand.grid(
  social_ties = 0:3,
  elec_comp   = c(0, 0.75),
  region = levels(droplevels(bwgdata$region))[1],
  purge       = median(bwgdata$purge, na.rm = TRUE),
  gwf_military= median(bwgdata$gwf_military, na.rm = TRUE),
  log_gdp     = mean(bwgdata$log_gdp, na.rm = TRUE),
  log_pop     = mean(bwgdata$log_pop, na.rm = TRUE),
  mid         = median(bwgdata$mid, na.rm = TRUE),
  civil_war   = median(bwgdata$civil_war, na.rm = TRUE),
  coups_cumulative = median(bwgdata$coups_cumulative, na.rm = TRUE),
  log_milex   = mean(bwgdata$log_milex, na.rm = TRUE),
  hhi = mean(bwgdata$hhi, na.rm=TRUE),
  log_milper  = mean(bwgdata$log_milper, na.rm = TRUE),
  media_diverse = mean(bwgdata$media_diverse, na.rm=TRUE),
  KEEP.OUT.ATTRS = FALSE, stringsAsFactors = FALSE
)

# Region Fixed Effect #

pred_data$region <- factor(pred_data$region, levels = levels(bwgdata$region))


# Pulling Count and Zero Coefficients #

all_names <- colnames(robust_vcov2)
count_names_full <- all_names[grepl("^count_", all_names)]
zero_names_full  <- all_names[grepl("^zero_",  all_names)]

count_names <- sub("^count_", "", count_names_full)
zero_names  <- sub("^zero_",  "", zero_names_full)

# Model Matrix #

X_count <- model.matrix(
  ~ social_ties * elec_comp + region + purge + gwf_military +
    log_gdp + log_pop + mid + civil_war + coups_cumulative +
    log_milex + log_milper + hhi,
  data = pred_data
)
Z_zero  <- model.matrix(~ mid + civil_war + media_diverse + log_milper, data = pred_data)

X_count <- X_count[, count_names, drop = FALSE]
Z_zero  <- Z_zero[,  zero_names,  drop = FALSE]

beta  <- coef(model_zinb2, model = "count")
gamma <- coef(model_zinb2, model = "zero")

# Computing Predictions #

eta_count <- as.numeric(X_count %*% beta)        # log mu
mu        <- exp(eta_count)                      # mu
eta_zero  <- as.numeric(Z_zero %*% gamma)        # logit(pi)
pi0       <- plogis(eta_zero)                    # pi
Ehat      <- (1 - pi0) * mu                      # expected counts

# Adding Robust Errors #

logtheta_name <- all_names[grepl("^log\\(theta\\)", all_names)]
p_all <- length(all_names)

grad_list <- vector("list", nrow(pred_data))
for (i in seq_len(nrow(pred_data))) {
  g_count <- Ehat[i] * X_count[i, ]
  g_zero  <- -mu[i] * pi0[i] * (1 - pi0[i]) * Z_zero[i, ]
  g <- numeric(p_all)
  g[match(count_names_full, all_names)] <- g_count
  g[match(zero_names_full,  all_names)] <- g_zero
  grad_list[[i]] <- g
}

vars <- vapply(grad_list, function(g) as.numeric(t(g) %*% robust_vcov2 %*% g), numeric(1))
se   <- sqrt(pmax(vars, 0))

# CIs #

lwr <- pmax(Ehat - 1.96 * se, 0)
upr <- Ehat + 1.96 * se

pred_plot <- pred_data %>%
  mutate(
    fit = Ehat,
    lwr = lwr,
    upr = upr,
    elec_comp = factor(elec_comp, levels = c(0, 0.75), labels = c("0", "0.75"))
  )

#################
### Figure 2 ####
#################

pd <- position_dodge(width = 0.3)

ggplot(pred_plot, aes(x = social_ties, y = fit, color = elec_comp, group = elec_comp)) +
  geom_point(size = 3, position = pd) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.15, position = pd) +
  labs(x = "Social Ties", y = "",
       color = "Electoral Competition") +
  theme_minimal(base_size = 12)


#########################
### Table 3 - Model 1 ###
#########################

model_praise <- glm(
  praise_binary ~ social_ties * elec_comp + region + purge + gwf_military +
    log_gdp + log_pop + mid + civil_war + coups_cumulative +
    log_milex + log_milper,
  data = bwgdata,
  family = binomial(link = "logit")
)

robust_vcov <- sandwich::vcovHC(model_praise, type = "HC0")
coeftest(model_praise, vcov = robust_vcov)

# Prediction #

pred_data <- expand.grid(
  social_ties = 0:3,
  elec_comp   = c(0, 0.75),
  region      = levels(bwgdata$region)[1],
  purge       = median(bwgdata$purge, na.rm = TRUE),
  gwf_military= median(bwgdata$gwf_military, na.rm = TRUE),
  log_gdp     = mean(bwgdata$log_gdp, na.rm = TRUE),
  log_pop     = mean(bwgdata$log_pop, na.rm = TRUE),
  mid         = median(bwgdata$mid, na.rm = TRUE),
  civil_war   = median(bwgdata$civil_war, na.rm = TRUE),
  coups_cumulative = median(bwgdata$coups_cumulative, na.rm = TRUE),
  log_milex   = mean(bwgdata$log_milex, na.rm = TRUE),
  log_milper  = mean(bwgdata$log_milper, na.rm = TRUE)
)

# Region Fixed Effect #

pred_data$region <- factor(pred_data$region, levels = levels(bwgdata$region))

# Model Matrix #

X <- model.matrix(
  ~ social_ties * elec_comp + region + purge + gwf_military +
    log_gdp + log_pop + mid + civil_war + coups_cumulative +
    log_milex + log_milper,
  data = pred_data
)

# Computing Predictions #

coefs <- coef(model_praise)
linpred <- X %*% coefs
se <- sqrt(diag(X %*% robust_vcov %*% t(X)))

pred <- plogis(linpred)
lwr  <- plogis(linpred - 1.96*se)
upr  <- plogis(linpred + 1.96*se)

pred_plot <- pred_data %>%
  mutate(
    fit = pred,
    lwr = lwr,
    upr = upr,
    elec_comp = factor(elec_comp)
  )

#################
### Figure 5 ####
#################

pd <- position_dodge(width = 0.3)
ggplot(pred_plot, aes(x = social_ties, y = fit, color = elec_comp, group = elec_comp)) +
  geom_point(size = 3, position = pd) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.15, position = pd) +
  labs(x = "Social Ties", y = "",
       color = "Electoral Competition") +
  theme_minimal(base_size = 12)

#########################
### Table 2 - Model 3 ###
#########################

model_zinb3 <- zeroinfl(
  mps_count ~ social_ties * execnat + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi| mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov3 <- sandwich::vcovOPG(model_zinb3)
coeftest(model_zinb3, vcov = robust_vcov3)

# Prediction #

pred_data <- expand.grid(
  social_ties = 0:3,
  execnat    = c(0, 1),
  region      = levels(bwgdata$region)[1],
  purge       = median(bwgdata$purge, na.rm = TRUE),
  gwf_military= median(bwgdata$gwf_military, na.rm = TRUE),
  log_gdp     = mean(bwgdata$log_gdp, na.rm = TRUE),
  log_pop     = mean(bwgdata$log_pop, na.rm = TRUE),
  mid         = median(bwgdata$mid, na.rm = TRUE),
  civil_war   = median(bwgdata$civil_war, na.rm = TRUE),
  coups_cumulative = median(bwgdata$coups_cumulative, na.rm = TRUE),
  log_milex   = mean(bwgdata$log_milex, na.rm = TRUE),
  log_milper  = mean(bwgdata$log_milper, na.rm = TRUE),
  hhi         = mean(bwgdata$hhi, na.rm = TRUE),
  media_diverse = mean(bwgdata$media_diverse, na.rm=TRUE),
  KEEP.OUT.ATTRS = FALSE, stringsAsFactors = FALSE
)

# Region Fixed Effect #

pred_data$region <- factor(pred_data$region, levels = levels(bwgdata$region))

# Pulling Count and Zero Coefficients #

all_names <- colnames(robust_vcov3)
count_names_full <- all_names[grepl("^count_", all_names)]
zero_names_full  <- all_names[grepl("^zero_",  all_names)]

count_names <- sub("^count_", "", count_names_full)
zero_names  <- sub("^zero_",  "", zero_names_full)

# Model Matrix #

X_count <- model.matrix(
  ~ social_ties * execnat + region + purge + gwf_military +
    log_gdp + log_pop + mid + civil_war + coups_cumulative +
    log_milex + log_milper + hhi,
  data = pred_data
)
Z_zero <- model.matrix(~ mid + civil_war + media_diverse + log_milper, data = pred_data)

X_count <- X_count[, count_names, drop = FALSE]
Z_zero  <- Z_zero[,  zero_names,  drop = FALSE]

# Computing Predictions #

beta  <- coef(model_zinb3, model = "count")
gamma <- coef(model_zinb3, model = "zero")

eta_count <- as.numeric(X_count %*% beta) 
mu        <- exp(eta_count)
eta_zero  <- as.numeric(Z_zero %*% gamma)
pi0       <- plogis(eta_zero)
Ehat      <- (1 - pi0) * mu

# CIs #

logtheta_name <- all_names[grepl("^log\\(theta\\)", all_names)]
p_all <- length(all_names)

grad_list <- vector("list", nrow(pred_data))
for (i in seq_len(nrow(pred_data))) {
  g_count <- Ehat[i] * X_count[i, ]
  g_zero  <- -mu[i] * pi0[i] * (1 - pi0[i]) * Z_zero[i, ]
  
  g <- numeric(p_all)
  g[match(count_names_full, all_names)] <- g_count
  g[match(zero_names_full,  all_names)] <- g_zero
  
  grad_list[[i]] <- g
}

vars <- vapply(grad_list, function(g) as.numeric(t(g) %*% robust_vcov3 %*% g), numeric(1))
se   <- sqrt(pmax(vars, 0))

lwr <- pmax(Ehat - 1.96 * se, 0)
upr <- Ehat + 1.96 * se

pred_plot <- pred_data %>%
  mutate(
    fit = Ehat,
    lwr = lwr,
    upr = upr,
    execnat = factor(execnat, labels = c("0","1"))
  )

#################
### Figure 3 ####
#################

pd <- position_dodge(width = 0.3)
ggplot(pred_plot, aes(x = social_ties, y = fit, color = execnat, group = execnat)) +
  geom_point(size = 3, position = pd) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.15, position = pd) +
  labs(x = "Social Ties", y = "",
       color = "HOG Nationalist") +
  theme_minimal(base_size = 12)

#########################
### Table 3 - Model 2 ###
#########################

model_criticism <- glm(
  criticism_binary ~ social_ties * execnat + region + purge + gwf_military +
    log_gdp + log_pop + mid + civil_war + coups_cumulative +
    log_milex + log_milper,
  data = bwgdata,
  family = binomial(link = "logit")
)

robust_vcov2 <- vcovHC(model_criticism, type = "HC0")
coeftest(model_criticism, vcov = robust_vcov2)

# Prediction #

pred_data <- expand.grid(
  social_ties = 0:3,
  execnat    = c(0, 1),
  region      = levels(bwgdata$region)[1],
  purge       = mean(bwgdata$purge, na.rm = TRUE),
  gwf_military= mean(bwgdata$gwf_military, na.rm = TRUE),
  log_gdp     = mean(bwgdata$log_gdp, na.rm = TRUE),
  log_pop     = mean(bwgdata$log_pop, na.rm = TRUE),
  mid         = mean(bwgdata$mid, na.rm = TRUE),
  civil_war   = mean(bwgdata$civil_war, na.rm = TRUE),
  coups_cumulative = mean(bwgdata$coups_cumulative, na.rm = TRUE),
  log_milex   = mean(bwgdata$log_milex, na.rm = TRUE),
  log_milper  = mean(bwgdata$log_milper, na.rm = TRUE)
)

# Region Fixed Effect #

pred_data$region <- factor(pred_data$region, levels = levels(bwgdata$region))

# Model Matrix #

X <- model.matrix(
  ~ social_ties * execnat + region + purge + gwf_military +
    log_gdp + log_pop + mid + civil_war + coups_cumulative +
    log_milex + log_milper,
  data = pred_data
)

# Computing Predictions #

coefs <- coef(model_criticism)
stopifnot(identical(colnames(X), names(coefs)))

linpred <- X %*% coefs
se <- sqrt(diag(X %*% robust_vcov2 %*% t(X)))

pred <- plogis(linpred)
lwr  <- plogis(linpred - 1.96*se)
upr  <- plogis(linpred + 1.96*se)

pred_plot <- pred_data %>%
  mutate(
    fit = pred,
    lwr = lwr,
    upr = upr,
    execnat = factor(execnat, labels = c("0","1"))
  )

#################
### Figure 4 ####
#################

pd <- position_dodge(width = 0.3)
ggplot(pred_plot, aes(x = social_ties, y = fit, color = execnat, group = execnat)) +
  geom_point(size = 3, position = pd) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.15, position = pd) +
  labs(x = "Social Ties", y = "",
       color = "HOG Nationalist") +
  theme_minimal(base_size = 12)

###########################
##### Online Appendix #####
###########################

### Table 4 ###
###############

model_beta <- betareg(
  trust_pct ~ social_ties + region + purge + gwf_military +
    log_gdp + log_pop + mid + civil_war + coups_cumulative +
    log_milex + log_milper,
  data = bwgdata
)

robust_vcov <- sandwich::vcovHAC(model_beta)
coeftest(model_beta, vcov = robust_vcov)

### Table 5 ###
###############

# Model 1 #

model_zinb <- zeroinfl(
  mps_count ~ social_ties * elec_comp + lag(mps_count) + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi| mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb)
coeftest(model_zinb, vcov = robust_vcov)

# Model 2 #

model_zinb2 <- zeroinfl(
  mps_count ~ social_ties * execnat + lag(mps_count) + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi| mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb2)
coeftest(model_zinb2, vcov = robust_vcov)

### Table 6 ###
###############

# Model 1 #

model_zinb <- zeroinfl(
  mps_count ~ com_mil_serv * elec_comp + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi| mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb)
coeftest(model_zinb, vcov = robust_vcov)

# Model 2 #

model_zinb <- zeroinfl(
  mps_count ~ mil_police * elec_comp + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi| mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb)
coeftest(model_zinb, vcov = robust_vcov)

# Model 3 #

model_zinb <- zeroinfl(
  mps_count ~ mil_eco_dummy * elec_comp + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi| mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb)
coeftest(model_zinb, vcov = robust_vcov)

### Table 7 ###
###############

# Model 1 #

model_zinb2 <- zeroinfl(
  mps_count ~ com_mil_serv * execnat + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi| mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb2)
coeftest(model_zinb2, vcov = robust_vcov)

# Model 2 #

model_zinb2 <- zeroinfl(
  mps_count ~ mil_police * execnat + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi| mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb2)
coeftest(model_zinb2, vcov = robust_vcov)

# Model 3 #

model_zinb2 <- zeroinfl(
  mps_count ~ mil_eco_dummy * execnat + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi| mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb2)
coeftest(model_zinb2, vcov = robust_vcov)

### Table 8 ###
###############

model_zinb <- zeroinfl(
  mps_count ~ social_ties * polity + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi| mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb)
coeftest(model_zinb, vcov = robust_vcov)

### Table 9 ###
###############

# ZINB Model #

model_zinb <- zeroinfl(
  mps_count ~ social_ties * elec_comp + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi | mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov_zinb <- sandwich::vcovOPG(model_zinb)
coeftest(model_zinb, vcov = robust_vcov_zinb)

# Negative Binomial Model #

model_nb <- glm.nb(
  mps_count ~ social_ties * elec_comp + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi,
  data = bwgdata
)

robust_vcov_nb <- sandwich::vcovHC(model_nb, type = "HC0")
coeftest(model_nb, vcov = robust_vcov_nb)

# AIC Comparison #

AIC(model_zinb, model_nb)

# ZINB Model 2 #

model_zinb2 <- zeroinfl(
  mps_count ~ social_ties * execnat + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi | mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb2)
coeftest(model_zinb2, vcov = robust_vcov)

# Negative Binomial Model 2 #

model_nb2 <- glm.nb(
  mps_count ~ social_ties * execnat + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi,
  data = bwgdata
)

robust_vcov_nb2 <- sandwich::vcovHC(model_nb2, type = "HC0")
coeftest(model_nb2, vcov = robust_vcov_nb2)

# AIC Comparison 2 #

AIC(model_zinb2, model_nb2)


### Table 10 ###
################

# Model 1 #

model_zinb <- zeroinfl(
  mps_count ~ social_ties * elec_comp + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi + media_censor | mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov_zinb <- sandwich::vcovOPG(model_zinb)
coeftest(model_zinb, vcov = robust_vcov_zinb)

# Model 2 #

model_zinb2 <- zeroinfl(
  mps_count ~ social_ties * execnat + region + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi + media_censor | mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb2)
coeftest(model_zinb2, vcov = robust_vcov)

### Table 11 ###
################

# Model 1 #

model_zinb <- zeroinfl(
  mps_count ~ social_ties + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi | mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov_zinb <- sandwich::vcovOPG(model_zinb)
coeftest(model_zinb, vcov = robust_vcov_zinb)

# Model 2 #

model_zinb2 <- zeroinfl(
  mps_count ~ social_ties * elec_comp + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi | mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov_zinb <- sandwich::vcovOPG(model_zinb2)
coeftest(model_zinb2, vcov = robust_vcov_zinb)

# Model 3 #

model_zinb3 <- zeroinfl(
  mps_count ~ social_ties * execnat + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi | mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb3)
coeftest(model_zinb3, vcov = robust_vcov)

### Table 12 ###
################

# Model 1 #

model_praise <- glm(
  praise_binary ~ social_ties * elec_comp + purge + gwf_military +
    log_gdp + log_pop + mid + civil_war + coups_cumulative +
    log_milex + log_milper,
  data = bwgdata,
  family = binomial(link = "logit")
)

robust_vcov <- sandwich::vcovHC(model_praise, type = "HC0")
coeftest(model_praise, vcov = robust_vcov)

# Model 2 #

model_criticism <- glm(
  criticism_binary ~ social_ties * execnat + purge + gwf_military +
    log_gdp + log_pop + mid + civil_war + coups_cumulative +
    log_milex + log_milper,
  data = bwgdata,
  family = binomial(link = "logit")
)

robust_vcov <- sandwich::vcovHC(model_criticism, type = "HC0")
coeftest(model_criticism, vcov = robust_vcov)

### Table 13 ###
################

model_zinb <- zeroinfl(
  mps_count ~ social_ties * execnat + region + anti_regime_mil + purge + gwf_military + 
    log_gdp + log_pop + mid + civil_war + coups_cumulative + 
    log_milex + log_milper + hhi | mid + civil_war + media_diverse + log_milper,
  data = bwgdata,
  dist = "negbin"
)

robust_vcov <- sandwich::vcovOPG(model_zinb)
coeftest(model_zinb, vcov = robust_vcov)

### Table 14 ###
################

model_criticism <- glm(
  criticism_binary ~ social_ties * execnat + region + anti_regime_mil + purge + gwf_military +
    log_gdp + log_pop + mid + civil_war + coups_cumulative +
    log_milex + log_milper,
  data = bwgdata,
  family = binomial(link = "logit")
)

robust_vcov <- sandwich::vcovHC(model_criticism, type = "HC0")
coeftest(model_criticism, vcov = robust_vcov)

### Table 17 ###
################

# Model 1 #

model_criticism <- glm(
  criticism_dummy ~ social_ties * execnat + region + purge + gwf_military +
    log_gdp + log_pop + mid + civil_war + coups_cumulative +
    log_milex + log_milper,
  data = bwgdata,
  family = binomial(link = "logit")
)

robust_vcov2 <- vcovHC(model_criticism, type = "HC0")
coeftest(model_criticism, vcov = robust_vcov2)

# Model 2 #

model_praise <- glm(
  praise_dummy ~ social_ties * elec_comp + region + purge + gwf_military +
    log_gdp + log_pop + mid + civil_war + coups_cumulative +
    log_milex + log_milper,
  data = bwgdata,
  family = binomial(link = "logit")
)

robust_vcov2 <- vcovHC(model_praise, type = "HC0")
coeftest(model_praise, vcov = robust_vcov2)
